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Abstract: Satellite observations of surface reflected solar radiation contain information 
about variability in the absorption of solar radiation by vegetation. Understanding the 
causes of variability is important for models that use these data to drive land surface fluxes 
or for benchmarking prognostic vegetation models. Here we evaluated the interannual 
variability in the new 30.5-year long global satellite-derived surface reflectance index data, 
Global Inventory Modeling and Mapping Studies normalized difference vegetation index 
(GIMMS NDVI3g). Pearson’s correlation and multiple linear stepwise regression analyses 
were applied to quantity the NDVI interannual variability driven by climate anomalies, and 
to evaluate the effects of potential interference (snow, aerosols and clouds) on the NDVI 
signal. We found ecologically plausible strong controls on NDVI variability by antecedent 
precipitation and current monthly temperature with distinct spatial patterns. Precipitation 
correlations were strongest for temperate to tropical water limited herbaceous systems 
where in some regions and seasons > 40% of the NDVI variance could be explained by 
precipitation anomalies. Temperature correlations were strongest in northern mid- to 
high-latitudes in the spring and early summer where up to 70% of the NDVI variance was 
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explained by temperature anomalies. We find that, in western and central North America, 
winter-spring precipitation determines early summer growth while more recent 
precipitation controls NDVI variability in late summer. In contrast, current or prior wet 
season precipitation anomalies were correlated with all months of NDVI in sub-tropical 
herbaceous vegetation. Snow, aerosols and clouds as well as unexplained phenomena still 
account for part of the NDVI variance despite corrections. Nevertheless, this study 
demonstrates that GIMMS NDVI3g represents real responses of vegetation to climate 
variability that are useful for global models. 

Keywords: GIMMS NDVI3g; climate-driven interannual variability; interference 


1. Introduction 

Over the past 30 years atmospheric CO 2 levels have been rising as a result of fossil fuel emissions. 
Superimposed on the atmospheric CO 2 growing trend are large interannual excursions (nearly 100% of 
the trend), and evidence strongly suggests that this interannual variability in the atmospheric CO 2 
growth rate is driven largely by the terrestrial biosphere [1,2]. Global carbon cycle models are 
challenged to reproduce this CO 2 growth rate variability, and explanations largely involve independent 
responses of terrestrial net primary production (NPP), heterotrophic respiration (RH), and fire 
emissions to climate anomalies such as those associated with ENSO events [3-6]. Understanding the 
potential for the terrestrial biosphere to mitigate or perhaps aggravate CO 2 accumulation in the 
atmosphere may be as important for predicting future climate change as projections of future fossil 
fuel emissions. Improved forecasts of future trends in atmospheric CO 2 will be possible, through the 
development of models that can realistically represent the processes controlling current and past 
terrestrial carbon fluxes [7]. 

Precise modeling of NPP is key in terrestrial carbon cycle models, because NPP not only represents 
the net uptake of carbon from the atmosphere by vegetation but also affects other carbon fluxes by 
providing the substrate for RH and fuel for combustion by fire. Terrestrial NPP has been showed to 
play a central role in determining the local and global CO 2 content of the atmosphere at temporal 
scales spanning hours [8] to epochs [9]. 

Long-term satellite observations of vegetation “greenness” provide a record of variability that 
greatly improves modeling NPP in terrestrial carbon cycle models. Satellite derived greenness refers to 
various surface reflectance indices that can be measured with satellite instruments. These indices 
express the amount of green vegetation absorbing solar radiation, a state that is often nearly 
proportional to photosynthetic productivity [10-12] and net ecosystem carbon exchange [13]. Global 
carbon cycle models often rely on satellite observations of vegetation greenness to either prescribe 
phenological variability [14] or to evaluate vegetation dynamics in prognostic models [15]. 

Polar orbiting satellites have provided global measurements of the greenness index Normalized 
Difference Vegetation Index (NDVI) at sub-monthly time steps and/or moderate spatial resolutions 
(>30 m) for decades. NDVI is commonly used to derive canopy solar radiation absorptance for 
estimating NPP [3,16,17]. It has also been used in global climate models to prescribe vegetation 
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influences on the hydrological and surface energy cycle [18-21]. NDVI has been shown to indicate 
variability in vegetation productivity in response to climate variability [22] and human management [23]. 
NDVI may not capture all the climate driven variability in NPP, however, and includes noise and biases 
that interfere with diagnosing true vegetation responses to climate [24]. Even different processing of the 
same NDVI data set can result in large differences in modeled global carbon fluxes [25]. 

The purpose of this study is to quantity how much of the interannual variability in a new NDVI data 
set reflects real responses of vegetation to interannual climate variability, and to evaluate whether the 
responses of NDVI to climate variability are meaningful enough to be of value for modeling 
vegetation activity. This NDVI data set is the latest version of Global Inventory Modeling and 
Mapping Studies (GIMMS) NDVI3g. It is derived from NOAA’s Advanced Very High Resolution 
Radiometer (AVHRR) satellite record, and is the longest global sub-monthly time series of greenness 
index available (July 1981-present), more than twice as long as those available from newer sensors 
such as NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS, February 2000-present). 
Since this NDVI data set is used as input to or for benchmarking global carbon, water and energy cycle 
models, it is important that we understand how much variability in vegetation growth is captured by 
this NDVI data set and if the causes of variability are linked to climate (precipitation, temperature, and 
solar radiation), disturbance (fires and large area outbreaks of pests), human management (e.g., 
irrigation and fertilization), or residual errors. Other studies have been published that reveal 
correlations between climate and NDVI [26-30], but none have combined analyses on monthly 
anomalies, lead time dynamics, cumulative climate effects and non-climate signal interference at 
global scales [31,32], and none have used global NDVI time series longer than 20 years. 

We explored the sources of variability in the GIMMS NDVI3g, in particular the ability of the 
GIMMS NDVI3g to capture vegetation responses to climate variability by: (1) evaluating the performance 
of GIMMS NDVBg anomalies in comparison to the NDVI anomalies from the more advanced MODIS 
instrument; (2) performing correlation analysis between the anomalies of GIMMS NDVBg and climate 
variables (precipitation and temperature) and comparing our results to previous studies; and (3) applying 
multiple linear regression to estimate the amount of GIMMS NDVBg variance that is driven by climate 
versus that caused by interference in surface reflectance signal (e.g., snow, aerosols, clouds, and residual 
effects). Our results corroborate and expand in detail previously reported responses of NDVI to climate 
variability. We argue that a significant part of the variability in GIMMS NDVBg reflects a true response 
of vegetation to climate variability at monthly to interannual time scales. Our analyses identified unique 
regional and seasonally specific dynamics in the response of vegetation to antecedent precipitation. 
Identification and quantification of the climate responses of vegetation represented in the GIMMS 
NDVBg data will lead to improved understanding and prediction of interannual variability of, among 
other things, terrestrial carbon fluxes and atmospheric CO 2 growth rate. 

2. Data Sets and Methods 

2.1. Data Sets 

The GIMMS NDVBg data set (here after GIMMS NDVI or NDVI unless specifically defined 
otherwise) was processed in a way consistent with and quantitatively comparable to NDVI derived 
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from improved sensors such as MODIS and SPOT-4 Vegetation, and was corrected for dropped scan 
lines, navigation errors, data drop outs, edge-of-orbit composite discontinuities, and other 
artifacts [33]. The data processing algorithm also minimized solar zenith angle effects introduced by 
orbital drift and effects of changes in the sun-target-sensor geometry, corrected for bias using 
Sea-viewing Wide Field-of-view Sensor (SeaWiFS) data, and corrected for volcanic stratospheric 
aerosol effects from the El Chichon (1982-1984) and Mt Pinatubo (1991-1993) volcanic 
eruptions [33,34]. Calibration and correction for orbital drift artifacts were accomplished by empirical 
mode decomposition (EMD); cloud contamination was fdtered when the NDVI values drop beyond 
95% confidence interval; and maximum value compositing (MVC) was used to minimize atmospheric 
(e.g., aerosols) and radiative geometry effects [33,34]. For MODIS NDVI, in addition to MVC a more 
sophisticated approach was used for atmospheric corrections: NDVI was derived from surface 
reflectances corrected for atmospheric effects [35]. The atmospheric correction algorithm of MODIS 
surface reflectances used an atmospheric radiative transport model and accounted for the effects of 
gaseous (03, O 2 , CO 2 , and water vapor) and aerosol scattering and absorption, adjacency effects 
caused by variation of land cover, surface/atmosphere Bidirectional Reflectance Distribution Function 
(BRDF), atmosphere coupling effects, and contamination by thin cirrus [36,37]. 

We acquired GIMMS NDVI at the native resolution of 0.083°, bi-monthly. Since long term 
observational global climate data sets are at much coarser spatial scales, we focused our analysis on 
spatial scales of 0.25° to 1°. We aggregated the data by averaging to 0.25°, 0.5°, and 1° spatial 
resolutions, and by averaging bi-monthly to produce monthly time steps, depending on the native 
resolutions of the comparison data sets. The temporal and spatial resolutions considered here reflect 
those typical to global carbon, water, and energy cycle models especially those coupled to atmospheric 
transport. Monthly MODIS Terra (MOD13C2 V5, 2001-2010) and Aqua (MYD13C2 V5, 2003-2010) 
NDVI were also acquired (http://reverb.echo.nasa.gov) and regridded by averaging from 0.05° to 0.25° 
for comparisons. 

We used the newly released monthly Global Precipitation Climatology Project (GPCP) V2.2 
(http://precip.gsfc.nasa.gov). This data set has complete global coverage from January 1979 to 
December 2010, and a native spatial resolution of 2.5° [38]. We disaggregated the 2.5° data to 0.5° 
with no interpolation, and then averaged 0.5° to obtain the 1° data. We assume that the scales of 
significant regional climate anomalies are on the order of or larger than 2.5° x 2.5°. For comparison 
we also obtained Tropical Rainfall Measuring Mission (TRMM, 1998-2010) 0.25° daily precipitation 
product (3B42 V6, http://mirador.gsfc.nasa.gov) with nominal spatial domain of 50°N-50°S, and 
aggregated it from daily to bi-monthly and monthly temporal resolutions. 

We constructed the temperature data set from the 0.5° Climate Research Unit (CRU CL 1.0, 
http://www.ipcc-data.org/obs/get_30yr_means.html) mean monthly climatology [39] and 2° Goddard 
Institute for Space Studies (GISS, 1200 km smoothing radius, http://www.giss.nasa.gov) monthly 
surface anomalies which extend for more than 130 years from 1880 to present [40]. The base period of 
the CRU climatology and GISS anomalies was 1961-1990 [39]. The derived temperature data set has 
a spatial resolution of 0.5° and a monthly temporal resolution. 

Monthly MODIS snow, aerosol and cloud data for 2000-2010 were obtained to evaluate how these 
variables interfere with estimates of vegetation variability derived from NDVI. The snow product is 
MOD10CM V5 (http://reverb.echo.nasa.gov) at 0.5° spatial resolution aggregated from 0.05° 
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resolution [41]. Aerosol and cloud data were regridded to 0.5° from the 1° MOD08_M3 V5.1 product 
(ftp://ladsweb.nascom.nasa.gov). All regridding from fine to coarser scales was accomplished by 
averaging finer data within the coarser grid cell. Regridding from coarse to finer scale was achieved by 
assigning the value of the coarser grid cell to the finer grid cells within. Interpolation was not used 
in regridding. 

2.2. Methods 

For all the time series of GIMMS and MODIS NDVI, climate variables and potential interfering 
elements, we generated their monthly anomalies by subtracting their respective 29-year (for GIMMS 
NDVI, GPCP precipitation and GISS temperature, 1982-2010), 13-year (for TRMM precipitation, 
1998-2010), 11-year (for snow, aerosols and clouds, 2000-2010), or 8-year (for MODIS Terra and 
Aqua NDVI, 2003-2010) averaged annual cycles following Equation (1) below: 

N 

7 monthly _ value(y, m) ^ 

monthly _ anomalv(y, m ) = monthly _ value{y, m ) - — 

N 

where y and m are the year and month of the monthly value, respectively, and N is the number of years 
of the data set. This procedure removes the seasonal cycles of the data sets and allows us to study the 
relationship between the interannual variabilities of NDVI and climate variables. The trends in the 
anomaly time series were not eliminated because they are considered part of the interannual 
variabilities, and they are strongly linked to climate change in some regions such as northern high 
latitudes [26]. The trends are generally very small (within ±0.01 yr -1 , [42-46]) compared to the 
interannual variabilities in NDVI, accounting for <10% of the anomalies. For temperature, the GISS 
anomalies were not directly used because their base period (1961-1990) is not consistent with that of 
the NDVI data (1982-2010). 

All the statistical analyses (correlations and regressions) were conducted on monthly anomalies 
unless specified (e.g., bi-monthly for GIMMS NDVI and TRMM precipitation anomaly correlations). 
NDVI and precipitation data used in the analyses were GIMMS NDVI3g and GPCP precipitation, 
respectively, unless specified. 

We conducted Pearson’s correlation analyses at the grid cell level between GIMMS NDVI and 
MODIS Aqua NDVI anomalies at 0.25° spatial resolution for 2003-2010, between GIMMS NDVI and 
GPCP precipitation anomalies at 1° spatial resolution for 1982-2010, and between GIMMS NDVI and 
GISS temperature anomalies at 0.5° spatial resolution for 1982-2010. Selected spatial resolutions 
reflect the compromise between fine spatial scale NDVI signal and the coarser scale climate data. 
Correlation between GIMMS and MODIS Aqua NDVI was conducted on anomalies of all months 
(N = 96 months of 8 years), and the results were compared with the correlation between MODIS Aqua 
and Terra NDVI anomalies (N = 96 months of 8 years). Aqua NDVI, instead of Terra NDVI, was used 
to correlate with GIMMS NDVI because of the sensor degradation issue detected on the Terra 
platform [47]. For GIMMS NDVI-GPCP precipitation and GIMMS NDVI-GISS temperature, 
correlation analyses were performed on the anomalies of both all months (N = 348 months of 29 years) 
and individual months (N = 29 years), each with varying climate lead times to examine how fast 
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vegetation reacts to climate anomalies. Only the correlation coefficients that were significant at 95% 
confidence level (p < 0.05, two-tailed) are reported. 

We recognized that precipitation can only approximate soil moisture limitations on vegetation. 
Actual soil moisture depends not only on precipitation but also on interception, evaporation, lateral 
and vertical runoff, and soil characteristic such as depth and water holding capacity. To account for an 
energy constraint on soil moisture we analyzed correlation patterns for individual months between the 
anomalies of GIMMS NDVI and a water availability index defined as GPCP precipitation minus 
potential evaporation derived from net radiation data (2000-2009, from CERES instruments, 
Wielicki et al. [48], http://ceres.larc.nasa.gov/) and GISS temperature following the Priestley-Taylor 
method [49]. The NDVI-water availability anomaly correlation patterns (not shown) were generally 
the same as those for GPCP precipitation alone though the NDVI-water availability anomaly 
correlations were degraded by the shorter time length of the data. 

We tested the assumption of Gaussian distribution of variance and possible bias caused by outliers 
using a non-parametric approach (Spearman’s rank correlation) on the NDVI-precipitation and 
NDVI-temperature correlations. We found good agreement between the results obtained from 
Pearson’s correlation and Spearman’s rank correlation (see Figure SI). All other results presented are 
based on the Pearson's correlation. 

Stepwise regression analysis was applied per grid cell to estimate the amount of NDVI variance that 
could be attributed to climate (antecedent precipitation and current monthly temperature), known 
interference for which data were available (current monthly precipitation, snow, aerosols and clouds), 
and a residual term that combined other unidentified sources of NDVI variability that might include 
disturbance, human management and underestimated errors. This analysis was done for each month to 
reveal the seasonality of the relationship of the variables. Only the coefficients of determination that 
were significant at 95% confidence level (p < 0.05) were reported. 

2.2. 1 . Correlation Analyses on the NDVI-Climate Anomalies of All Months 

The NDVI and climate anomalies were aggregated to a time series each (N = 348 months of 
29 years), and correlation was carried out for climate lead times up to 6 months for precipitation and 
up to 4 months for temperature. For the lead 0 case, correlation analysis was conducted on the full time 
series of NDVI and precipitation or temperature anomalies. For lead k (precipitation or temperature 
leading NDVI by k months), we dropped the first k values from the NDVI anomaly array and the last k 
values from the precipitation or temperature anomaly array and conducted the correlation analysis. 

We evaluated whether finer temporal and spatial resolutions influenced the strength or spatial 
distribution of correlations between NDVI and precipitation anomalies by comparing results from the 
1°, monthly NDVI and GPCP precipitation data which we used for the remainder of our analyses with 
results from 0.25°, bi-monthly NDVI and TRMM precipitation. The spatial distributions of 
correlations were largely unaffected by the temporal and spatial scale changes though NDVI-TRMM 
precipitation correlations were weaker due to the shorter time span of the TRMM data. 

Since we observed significant positive NDVI-precipitation correlations for multiple lead times, we 
further examined the correlations between NDVI anomaly and precipitation anomaly cumulated for 
varying lead time durations. We applied correlation analysis to NDVI anomaly of the current month 
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and precipitation anomaly summed over a duration ranging from 1 to 10 months (cumulative 
precipitation anomaly hereafter, N = 348 - duration + 1). Details about correlation analysis for varying 
durations can be found in Wang et al. [50]. Briefly, for duration of k months, the cumulative 
precipitation anomaly was the sum of precipitation anomalies over a k-month period corresponding to 
lead = 0 to k-1 months. The results of these analyses were presented in detail later but they showed 
negative correlations between NDVI anomalies and precipitation anomalies at 1 -month duration 
(/'. e. , lead 0) which were not ecologically meaningful and were likely the result of cloud interference 
and were thus not used in further analyses. The highest correlations were found between anomalies in 
prior 6 months cumulative precipitation and NDVI. From these results we identified 6-month duration 
from lead 1 to 6 months (lead 1 + 2 + , ... , + 6) to be the optimum duration to represent the total 
effects of precipitation on vegetation growth, and this 6-month duration was applied to correlation 
analyses on NDVI-cumulative precipitation correlations for individual months in Section 2.2.2. 

When dealing with correlations between two time series (e.g., the NDVI and precipitation 
anomalies of all months) we have to account for the effect of temporal autocorrelation (also called 
serial correlation), otherwise the significance of correlation would appear to be higher than it actually 
is (z'.e., giving a smaller p value). Ordinary correlation assumes that the observations of a variable are 
independent of each other. However, this assumption is generally invalid for time series. In a time 
series, very often the observation at time i is closely associated with the observations at time 
/- 1, i - 2, ... , meaning that temporal autocorrelation exists in the time series. For example, NDVI 
anomaly in June generally does not deviate much from NDVI anomaly in May. For a time series with 
N observations X; (/= 1, 2, ... ,N), if there is a significant correlation between the subsets 
X; (z'= 1, 2, ... , N—k) and Xi+k(i + k= \+k,2+k, , N), this correlation is called the Ath-order 

temporal autocorrelation and the associated correlation coefficient is the Ath-order temporal 
autocorrelation coefficient of this time series [51]. Few previous studies have corrected their 
correlation results for temporal autocorrelation. Without such correction, the strength of the 
correlations would be overestimated and the interpretation of the results would be misleading [51]. 

For our correlation analyses between NDVI and climate variable anomalies of all months, we 
accounted for the effect of the lst-order temporal autocorrelation following the procedure by Dawdy 
and Matalas [51]. Ideally higher orders of temporal autocorrelation should also be accounted for until 
the temporal autocorrelation of at least one of the two time series becomes insignificant. However, this 
process is very laborious given the large number of land grid cells. Practically only the lst-order 
temporal autocorrelation is considered [52]. This is reasonable because for NDVI and climate 
variables the correlation between values of current month and last month is likely the strongest [53]. 

Briefly, 3 steps [51] were applied to each land grid cell to estimate the correlation coefficient 
(rvalue) and significance of the correlation coefficient (p* value) that is corrected for the lst-order 
temporal autocorrelation. Take the GIMMS NDVI-GPCP precipitation anomaly correlation at lead 0 
as an example, so here the two time series are 1982-2010 monthly NDVI and precipitation anomalies 
(N = 348). First, the correlation coefficient (/-ndvi precipitation) and the associated p value were calculated. 
Second, the lst-order temporal autocorrelation coefficients z-ndvi(I) and r P reci P itation(l) of the NDVI 
anomaly time series and the precipitation anomaly time series, respectively, were computed. Finally, if 
both rNovi(l) and r P reci P itation(l) are significant at 95% confidence level, we derived the effective 
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number of degrees of freedom («*) and the adjusted p value (p*) from the sample size N and the 
correlation coefficients above using Equations (2) — (4): 


* _ „ v 1 fvDF/0) X r precipitation 0 ) 
1 + r NDVI (1) X r precipitation O') 

(2) 

1 n * 


^ NDVI-precipitation /N , L 2 

y ^ NDVI-precipitation 

(3) 

p* = 2 x [1 - \‘^f(t)dt] 

(4) 


here n = N - 2, and /(t) is the probability density function of Student’s /-distribution. If either /-ndvi(I) 
or r P recipitation(l) is insignificant at 95% confidence level, no correction is needed and p* is the same as 
p. The NDVI-precipitation anomaly correlation is considered significant if the p * value is below 0.05. 
In summary, the correlation coefficient rNovi precipitation remains the same but the significance of the 
correlation coefficient is degraded ( i.e.,p * > p) during this process. 

2.2.2. Correlation Analyses on the NDVI-Climate Anomalies of Individual Months 

Correlation analysis was also conducted for each month of the year (January to December) for 
climate lead times up to 7 months for precipitation and up to 4 months for temperature (N = 29). 
Correlations generally deteriorated after ~4 months, but in some places NDVI-precipitation 
correlations remained significant positive for lead times up to 6 months. Therefore, similar to the 
correlation of all-month data, for each month we accounted for the long-lasting effect of precipitation 
on NDVI by further analyzing the correlations between NDVI anomalies and precipitation anomalies 
cumulated over a 6-month duration from lead 1 to 6 months (N = 29), and the correlation results will 
be used in Section 2.2.3 to estimate the fraction of NDVI variance driven by climate. 

The NDVI-temperature correlations of northern high latitudes need to be interpreted with caution 
for the winter and spring months due to the interference of snow. In warmer winters and springs, there 
is less snow cover and the land surface produces a higher NDVI signal. However, this does not 
necessary mean a higher rate of vegetation growth. We detected potential snow interference in the 
NDVI-temperature correlations when NDVI-snow correlations were significantly negative (p < 0.05, 
N = 132 for all months and N = 1 1 for individual months). For both all months and individual months 
we corrected snow interference by setting the NDVI-temperature correlation of a grid cell to zero if 
significant negative NDVI-snow correlation was observed for this grid cell. 

2.2.3. NDVI Variance Explained by Climate 

The fraction of the NDVI variance explained by climate (R 2 ) can be approximately computed from 
the correlation coefficients (calculated in Section 2.2.2) of NDVI-cumulative precipitation (lead 1-6) 
and NDVI-temperature at lead 0 corrected for snow interference. NDVI-cumulative precipitation 
correlation was used to account for the long-lasting effect of precipitation on vegetation, and 
NDVI-temperature correlation at lead 0 was selected because it’s the strongest among all the lead 
times we tested (from 0 to 4 months). 
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For each land grid cell in each month, if only one of the two correlation coefficients (cumulative 
precipitation or temperature) was significant positive, R 2 was the square of the significant positive 
correlation coefficient; if both were significant positive, R 2 was essentially the coefficient of 
determination of the multiple linear regression of NDVI on cumulative precipitation and temperature 
and could be calculated using the following equation [54]: 


Y 

j^2 NDVI — precipitation 


Y — Z/ 

ND VI -temperature ND VI — precipitation 




1 — r. 


precipitation— temperature 


( 5 ) 


where r is the correlation coefficient. For all other cases, R 2 was zero. Only significant positive correlations 
were considered here because they were ecologically meaningful: increased vegetation activity with 
warmer temperatures in mid- and high- latitudes and with higher rainfall in water limited regions. 


2.2.4. NDVI Variance Explained by Accounting for Atmospheric Interference 


After estimating the variance in NDVI caused by climate variability we briefly addressed the 
variance not caused by climate. Snow, aerosols and clouds among others can interfere with detection 
of vegetation greenness by satellites. The GIMMS NDVI has been adjusted to account for many of 
these effects, but such effects have not been Hilly eliminated. The need for snow correction in the 
interpretation of the temperature correlations is one example of interfering elements that might 
contribute to NDVI variance. Precipitation of current month was considered as another potential 
interfering element because it was negatively correlated with NDVI in many regions (details will be 
presented in Section 3.2) which is not ecologically meaningful. 

To examine how much of the NDVI variance was explained by the effects of snow, aerosols, clouds 
and precipitation of the current month, forward stepwise multiple linear regression analysis [55] was 
applied to each grid cell and each of the 12 months. NDVI anomaly was the dependent variable, and 
anomalies of snow, aerosols, clouds and precipitation of the current month were independent variables 
(N = 1 1, 2000-2010). All of the independent variables included in the regression models were 
negatively correlated with NDVI anomalies. More than 85% of the final regression models (p < 0.05) 
were found to have only one independent variable, so the multicollinearity problem was negligible and 
no correction was needed. 


2.2.5. NDVI Variance Unexplained 

Even after taking into account climate and the interference identified above, there was still some of 
the NDVI variance that was not explained. To determine the extent of this fraction, forward stepwise 
multiple linear regression analysis was employed to the anomalies of NDVI, cumulative precipitation 
(lead 1-6), temperature, snow, aerosols, clouds and precipitation of the current month (lead 0) on all 
land grid cells where data were available, with NDVI anomaly as the dependent variable and the rest 
as independent variables (N= 11, 2000-2010). The coefficient of determination (R 2 ) of the regression 
represents the percentage of NDVI variance driven by climate anomalies and atmospheric interference 
listed above, and the fraction of the NDVI variance unexplained was estimated as 1 - R 2 . More than 
79% of the final regression models (p < 0.05) had only one independent variable, so again the 
multicollinearity problem was negligible. 
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3. Results and Discussion 

3.1. Comparisons with MODIS 

Though the GIMMS AVHRR data span almost three times the length of MODIS, the latter was 
designed to measure specific narrow reflectance bands, to have explicit atmospheric corrections, to be 
highly calibrated, and to minimize sun-surface-platfonn radiative geometry effects. Many design 
features of MODIS were implemented to overcome problems originally identified in the AVHRR 
record. The new version of the GIMMS NDVI data set represents the culmination of 3 decades of 
experience with AVHRR NDVI focused on identifying and minimizing known problems. To gain 
insight about the robustness of the GIMMS NDVI for analysis of its longer time period we compared 
GIMMS and MODIS monthly NDVI anomalies for the 2003-2010 Aqua period. Seventy six percent 
(76%) of land grid cells showed statistically significant (p < 0.05) positive correlations between 
GIMMS NDVI and MODIS Aqua NDVI anomalies (Figure S2a). This number was close to the 
percentage (87%) of land grid cells that had statistically significant (p < 0.05) positive correlations 
between the anomalies of MODIS Aqua and Terra NDVI (Figure S2b). Correlations between GIMMS 
and MODIS Aqua were generally high except over tropical forests, eastern temperate regions and a 
band across boreal Asia (Figure S2a). These same regions also showed lower correlations between 
Terra and Aqua NDVI (Figure S2b), indicating high uncertainty in the interannual NDVI signal in 
these regions for all sensors. Gallo et al. [56] and Fensholt and Proud [43] reported similar agreement 
between older versions of GIMMS NDVI and MODIS NDVI. 

These comparisons build confidence that at the scales of our analysis the GIMMS NDVI is 
comparable to the more advanced but shorter time length MODIS data. GIMMS NDVI should perform 
nearly as well as MODIS in detecting climate driven variability but in the case of GIMMS NDVI over 
the longer 30.5-year record. 

3.2. Precipitation Controls on NDVI Variability 

To characterize the response of vegetation measured as NDVI to precipitation we examined what 
regions and what latencies occur in the correlation between NDVI and precipitation. For each grid cell 
we aggregated all monthly NDVI anomalies (N = 348, January 1982 to December 2010) and identified 
significant correlations with precipitation anomalies at various precipitation anomaly lead times. We 
found that globally positive correlations were the strongest and most spatially extensive at lead 1 
(Figure la), and locally significant positive correlations could last up to a lead of 6 months. The 
patterns of positive correlations (Figure la) were closely associated with shrub, grass, and crop land 
cover types south of ~60°N (Figure lb, 2001 MODIS land cover, MCD12C1 V5, UMD Scheme, [57]), 
higher fractional herbaceous cover (Figure lc, MODIS Vegetation Continuous Fields MOD44B V5, [58]), 
and higher aridity index (Figure Id, [59]). Generally NDVI-precipitation anomaly correlations were 
strongest for grasslands, shrublands, croplands, and savannas and lowest for forests. The lower 
correlations between NDVI and precipitation in forests probably reflect in part lower water limitations 
compared to the regions of high correlations. 
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NDVI correlations with cumulative precipitation summed over various lead times were examined to 
account for the long lasting effects of precipitation. Figure 2a-c shows the spatial patterns in 
correlations with cumulative precipitation of all months (N = 348 - duration + 1) for durations of 1 to 
3 months (see Figure S3 for correlations for durations up to 10 months). The number of global land 
grid cells with significant positive correlations increased markedly from 1 -month duration (lead 0) to 
2-month duration (lead 0+1) and leveled off at 7-month duration (lead 0 + 1 + 6), and 80% of 

the land grid cells that have significant positive correlations at 7-month duration were significant at 
2-month duration (Figure 2d). The largest increase in the number of pixels with significant positive 
correlations occurred at lead 1 . Cumulative precipitation anomalies at 1 -month duration (lead 0) were 
negatively correlated with NDVI anomalies in some boreal regions and in part of Southeast Asia and 
South America (Figure 2a). Such negative correlations are not ecologically meaningful and thus 
precipitation of the current month was treated as an interfering element in the NDVI signal 
in Section 3.5. 

Results from lead correlation analyses for individual months showed that globally the pattern of 
strongest correlations at lead 1 was valid for all 12 months of the year, and the strongest correlations 
generally occurred during each region’s respective rainy season (Figure 3). For example, in central 
North America, central Eurasia, India and the Sahel, the strongest correlations with NDVI anomalies 
in June, July, August, September and October were for precipitation anomalies in May, June, July, 
August, and September, respectively (Figure 3f-j). Another example is eastern and southern Africa 
and Australia, where NDVI anomalies in each month of December to April were most strongly 
correlated with precipitation anomalies 1 month before (Figure 31 and 3a-d). Figure 4 shows the 
percentage of global land grid cells that have positive correlations (p < 0.05) between NDVI and 
precipitation anomalies as a function of lead times from 0 to 6 months for each month of the year. 
Often correlations were negative at lead 0 for some areas (Figure 2a), but strong positive correlations 
emerged at leads of 1 and 2 months and peaked at lead 1 with 8-12% of total land grid cells having 
significant positive NDVI-precipitation correlations for all 12 months of the year (Figure 4). Longer 
leads showed diminishing correlations. These global patterns are consistent with regional studies using 
either NDVI or enhanced vegetation index (EVI) and either ground- or satellite -based climate 
data [50,60-62], and with the global analysis of Lotsch et al. [32]. 

Despite the general global pattern of strongest correlations at 1 -month lead, in some temperate and 
subtropical herbaceous dominated regions NDVI responses to precipitation lead times varied seasonally, 
reflecting the dynamics of soil moisture availability. Prominent examples are described below. 

Correlation patterns over western and central North America exhibited seasonally varying 
dependence of NDVI on precipitation (Figures 5a and 3e-h): NDVI anomalies in May were most 
strongly correlated with winter precipitation anomalies; in June and July significant correlations were 
observed between NDVI and recent precipitation with some residual dependence of NDVI on 
winter/spring precipitation; by August only most recent precipitation history showed signification 
correlations with NDVI. These results are consistent with expected soil moisture storage dynamics: 
spring and early summer plant growth is dependent upon precipitation received during winter when 
precipitation exceeds evaporation and the excess precipitation is stored in snow and soils, whereas in 
late summer soil moisture from spring is largely depleted and thus plant growth depends on most 
recent precipitation. Although not pointed out in Mendez-Barroso et al. [61], similar patterns can also 
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be found in their results from the station-level time series of absolute EVI and precipitation values in 
northwestern Mexico which is among the regions presented here. 

Figure 4. Percentage of land grid cells that show significant {p < 0.05) positive correlation 
between GIMMS NDVI and GPCP precipitation anomalies for lead times ranging from 0 
to 6 months for 12 months of the year (N = 29, 1982 to 2010). For all 12 months, 1-month 
lead has the highest land coverage of significant positive NDVI-precipitation correlations. 
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An example of a different seasonal dynamic in NDVI responses to precipitation was seen in 
sub-tropical regions such as Australia. During the rainy season correlations were strongest for recent 
precipitation (October-May, Austral spring-fall), while in all dry season months (e.g., June-August, 
Austral winter) NDVI anomalies were correlated most strongly with previous rainy season 
precipitation (Figure 5b, Figure 3 and Figure S4). In other words, correlations with NDVI anomalies 
during the dry season were highest for precipitation in the previous February- April wet season. 
Another example of this pattern was found in Southern Africa (Figure S4). 
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3.3. Temperature Controls on NDVI Variability 

Combining all monthly NDVI and temperature anomalies (N = 348, January 1982 to December 2010) 
produced positive correlations between NDVI and temperature anomalies in northern mid- and 
high- latitudes (Figure 6). Zhou et al. [63] reported similar results in their northern latitude study. 
Correlations with temperature were highest at lead 0 (Figure 6). Central Europe shows the strongest 
temperature effect (Figure 6a) that lasts for 3 to 4 months (Figure 7), while in the other northern mid- and 
high- latitude regions the temperature effect lasts for only 1 to 2 months (Figure 7). This is likely a result of 
the much shorter snow cover duration in central Europe compared to other regions in the same or higher 
latitudes. MODIS Terra snow cover data show that in March snow cover has almost completely 
disappeared in central Europe while in other regions in the same or higher latitudes snow cover remains 
above 70%. The shorter snow cover duration in central Europe is also reflected in the earlier start of the 
growing season in central Europe compared to other regions in the same or higher latitudes (Figure 8). 

Seasonal positive correlations between NDVI and current temperature anomalies were evident and 
supported ecological interpretations. Positive correlations emerged in spring (March) and progressed 
northward into June (Figure 7), consistent with the field observations in northeastern Siberian tundra 
by Blok et al. [64] showing that early summer temperature is the most important factor influencing 
vegetation growth in mid- and high- latitudes. The spring-early-summer northward wave of 
correlation matches closely with the progress of the MODIS-derived start of the growing season 
(MOD12Q2, [65,66], Figure 8). A similar but weaker wave of correlation moved south from 
mid-summer to early fall (July-September, Figure 7). 

Figure 8. Start of growing season for cycle 1 in the Northern Hemisphere. Vegetation 
phenology data are the V4 MODIS Land Cover Dynamics (MOD12Q2, [65,66]) product from 
http://duckwater.bu.edu/lc/datasets.html. Unit: month in the year (Spatial resolution: 0.5°). 



3.4. Climate Driven Variance in NDVI 

As expected, we found that precipitation was rarely negatively correlated (p < 0.05) with NDVI for 
leads of 1 month or more. Temperature was always positively correlated with NDVI anomalies at a 
lead of 0 month except where temperature and precipitation anomalies were negatively correlated. The 
latter were water limited regions and low precipitation was associated with higher temperature and 
drought. Combining the variance explained by positive correlations with cumulative antecedent 
precipitation (1-6 months lead) and by positive correlations with temperature at lead 0 (excluding 
negatively correlated NDVI-snow grid cell-months) following the method described in Section 2.2.3 and 
Equation (5), we estimated the variance in monthly NDVI explained by climate (Figure 9). Percentage of 
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monthly NDVI variance ( R 2 ) driven by climate ranged from about 30% to 70% in some boreal regions of 
Eurasia and North America in spring and summer due to temperature, and in central and western North 
American, central Eurasia, Australia, Africa, and South America due to precipitation. 

Our approach to characterizing climate effects on NDVI was simple in that it considered only 
monthly precipitation and temperature (though we have looked at bi-monthly precipitation and NDVI 
and found little difference in results) and it used a single variant approach. Inclusion of the 
Priestley-Taylor energy constraint did not improve correlations with NDVI compared to 
precipitation-NDVI correlations. We have also applied a multi-variant approach but the results did not 
change our overall conclusions and were more difficult to interpret. 

Our simple statistical approach was applicable because the regions influenced by temperature 
versus precipitation tended to be distinct: temperature correlations were strongest above 30°N latitude 
while precipitation correlations were strongest for herbaceous vegetation, < 1000 mm per year 
precipitation below 60°N latitude. In addition we evaluated the correlations between temperature and 
precipitation (not shown) which showed little interactions in terms of effects on NDVI though the 
analysis helped to explain noise and bias in the NDVI. 

3.5. Influence of Snow, Aerosols, Clouds and Precipitation of the Current Month on NDVI Variability 

Some of the important conditions that interfere with satellite measurements of vegetation activity 
using NDVI include snow cover, atmospheric aerosols, and clouds. Precipitation of the current month 
showed negative correlations with NDVI anomalies in some regions (Figure 2a) and was also treated 
as an interfering element here. The readily available MODIS snow, aerosol and cloud products and 
precipitation data motivated exploration into how much variance in the NDVI signal is caused by these 
interfering elements that linger in the GIMMS NDVI3g. 

Figure 10 is an example month (September) from the series of monthly maps in Figure S5 that 
compare total monthly NDVI variance (Figures 10a and S5a), fractions of NDVI variance contributed 
by climate (Figures 10b and S5b, same as Figure 9g) and by atmospheric interference (Figures 10c 
and S5c), and fraction of unexplained NDVI variance (Figures lOd and S5d). Atmospheric interference 
includes snow, aerosols, clouds and precipitation of the current month all of which are negatively 
correlated with NDVI at lead 0 month (Section 2.2.4). Unexplained NDVI variance represents the 
residual NDVI variance that is not significantly correlated with cumulative precipitation (1-6 months 
lead) or current (lead 0) temperature, snow, aerosols, clouds, and precipitation anomalies (Section 
2.2.5). Generally NDVI variance was <0.01 (Figures 10a and S5a), which translates to a standard 
deviation of <0.1. In other words NDVI anomalies were less than 0.1 for most land grid cells most of 
time. Figure 10c suggests that interference by snow, aerosols, clouds and precipitation remains in the 
NDVI signal after the GIMMS NDVI3g processing. For some regions and some periods (e.g., North 
America, Europe and central Asia in winter months), more than 80% of the NDVI variance was 
associated with these atmospheric and surface interfering elements. Some notable patterns in the 
monthly maps in Figure S5 include (1) snow in April-May and climate (temperature) in June 
explained a large fraction of NDVI variance at high northern latitudes; and (2) the variability explained 
in the Amazon was mostly from clouds in February during the rainy season and from aerosols in 
September caused by fires (Figure 10c). 


Figure 9. Fraction of NDVI variance explained (R 2 ) by climate (cumulative precipitation for lead 1-6 and temperature at lead 0 corrected for 
snow interference, both positively correlated with NDVI with p < 0.05, 1982-2010) for the 12 months of the year (a-1) (Spatial resolution: 0.5°). 
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Figure 11. Percentage of vegetated grid cells in each Transcom region (a-k) with 
significant correlations {p < 0.05) between NDVI anomalies and anomalies in climate 
(cumulative precipitation of lead 1-6 and temperature at lead 0 corrected for snow 
interference, both positively correlated with NDVI, black solid line) and in current monthly 
precipitation, snow, aerosols and clouds (negatively correlated with NDVI with p < 0.05, 
excluding grid cells with significant positive NDVI-climate correlations, black dashed 
line), and percent remaining vegetated grid cells showing no correlations with the variables 
tested (unknown, gray dotted line). They sum up to 1 for each month in each region. 
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Even the combination of climate and the interference listed above cannot, folly explain all the NDVI 
variance. In regions such as Europe, southeastern China, parts of South America, and sometimes parts of 
boreal Asia, a significant fraction of NDVI variance remains unexplained (Figures lOd and S5d). 

Figure 1 1 summarizes regional (regions as defined in Transcom [67]) and seasonal patterns in % of 
land grid cells in which NDVI anomalies (1) are significantly correlated with climate (cumulative 
precipitation of lead 1-6 months and temperature at lead 0, all positive), (2) are accounted for variance 
(from snow, aerosols, clouds and precipitation of the current month combined, all negatively 
correlated with NDVI anomalies), and (3) show no correlations with the variables tested in this study. 
Significant climate correlations are generally most prevalent during the growing season (e.g., boreal 
and temperate regions, Southern Africa). Fraction of land grid cells in tropical regions (Tropical 
America and Asia) showing significant positive correlations with climate is low with little seasonal 
variability. Australia had high correlations with climate (precipitation) for all 12 months. Fraction of 
land grid cells with no significant climate or interference correlations was high in non-growing season 
in boreal and temperate regions. Tropical regions and temperate South America had high fractions of 
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land grid cells with no significant correlations with the variables tested. These results (Figures 9-1 1 
and S5) show when and where NDVI was more reliable at capturing vegetation responses to 
climate anomalies. 

3. 6. Uncertainties in the Data and Results 

All correlations and regressions presented here were those with probabilities >0.95 that were 
obtained solely by random chance. We tested the assumption of normality of distributions and the 
potential for bias caused by outliers by performing non-parametric Spearman’s rank correlation and 
found little influence on the results and conclusions. We tested for potential bias as a result of temporal 
and spatial aggregation by using a different precipitation data set (TRMM) with finer temporal and 
spatial resolutions (bi-monthly and 0.25° as opposed to monthly and 2.5° for original GPCP) and 
found that it did not change the results or conclusions. We tested to a limited extent the assumption 
that precipitation was an adequate surrogate for water limitation and found the assumption to be 
supported. We identified regions and months where correlations were most robust as well as 
where/when uncertainty was large. We evaluated the GIMMS NDVI against two other satellite data 
(MODIS Aqua and Terra) and identified where all three products were in agreement and where 
satellite NDVI data in general were not adequate to address vegetation responses to climate. 

There was still a significant fraction of NDVI variance unexplained in some regions and some 
periods. The unexplained NDVI variance may be partly due to the uncertainty in the NDVI and 
climate data sets. It could, however, be a result of the variability in other drivers of NDVI that we 
didn’t account for in this study (see the next section). 

3. 7 . Significant Drivers of NDVI Variance 

Other than climate driven variability, we have not taken into account processes that are likely to 
also influence NDVI variability, such as stimulation of plant growth caused by increasing atmospheric 
CO 2 , land use change such as deforestation and biomass burning, irrigation and fertilization, and N 
deposition. Particularly, land management practices such as irrigation and fertilization contribute to a 
significant fraction of the NDVI variability in croplands [68]. 

From our analysis we identify a number of significant causes of NDVI variance and their likely 
effects on NDVI. 

(1) Climate driven variability in NDVI as quantified in this study. Northern latitude NDVI 
anomalies were positively correlated with temperature. NDVI anomalies in temperate to 
tropical semi-arid and arid regions were positively correlated with precipitation with regionally 
varying seasonal dynamics. 

(2) Climate driven variance missed by our simple representation of climate variability (e.g., soil 
moisture, solar radiation, vapor pressure deficit (VPD), and evaporation) and by our linear 
spatial averaging. 

(3) Atmospheric interference from snow, aerosols and clouds reported here are negatively 
correlated with NDVI. Though these are addressed in part by the GIMMS NDVI3g processing, 
we found some residual negative correlations with these variables. 
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(4) Variance caused by sensor-sun-surface radiation geometry. This uncertainty is also addressed 
as part of the GIMMS NDVI3g processing prior to our use. 

(5) Other mechanistic processes not accounted for: iixigation, fertilization, fires and land degradation. 

(6) Atmospheric and surface interfering factors not accounted for in our analysis (e.g., water vapor 
and cloud type) and from aerosol, water vapor, and cloud effects not captured in the GIMMS 
NDVI3g processing. 

3.8. Future Work 

Our results support the use of satellite data for the representation of vegetation responses to climate 
variability in models. A few studies have confirmed that vegetation anomalies observed from satellites 
are associated with regional climate anomalies and are large enough to explain the anomalies in 
measured and modeled carbon and water fluxes [69,70], However, some global scale studies show 
little correspondence between satellite -based estimates of FPAR (fraction of photosynthetically active 
radiation) and anomalies in global carbon fluxes from atmospheric inversions [25,71]. Our next work 
will be to characterize the impacts of climate driven NDVI anomalies on modeled carbon fluxes and to 
evaluate these impacts using independent methods such as eddy covariance data and atmospheric CO 2 
inversions. We expect that variances of up to 0.01 in the NDVI signal in some regions during growing 
season can represent relative standard deviations of up to 10% in the FPAR which itself is proportional 
to NPP in some global carbon cycle models [17]. Preliminary evaluation of how this amount of climate 
driven variability influences carbon fluxes (not shown) suggests that climate driven variability in 
NDVI identified in this study could produce variability in simulated NPP of 20% or more in some 
regions and seasons. Examples of this level of variability in gross primary production and leaf area 
index have been reported at eddy covariance sites [72,73], and further work evaluating model 
responses to climate driven NDVI variability against other independent flux data is underway. 

4. Conclusions 

Our analyses confirm that GIMMS NDVI3g (Global Inventory Modeling and Mapping Studies 
normalized difference vegetation index, version 3) captures vegetation responses to climate nearly as 
well as more advanced sensors such as MODIS (Moderate Resolution Imaging Spectroradiometer) 
Aqua, in agreement with results from other studies [43,56]. The 30.5-year GIMMS AVHRR 
(Advanced Very High Resolution Radiometer) NDVI3g record thus provides a robust time series to 
explore causes of variations in land surface vegetation greenness as expressed by NDVI. 

In addition to taking advantage of the longer time series we also provided a more thorough 
month-by-month analysis of drivers of NDVI variability than previous global studies. Temperate and 
subtropical arid and semi-arid regions and northern high latitude spring and summer showed strong 
interannual climate driven variability in NDVI. Seasonal variability in the importance of antecedent 
precipitation represented seasonal variability in available soil moisture as in central North America 
where dependence of vegetation growth on water stored from the winter gradually diminishes into the 
summer. These patterns of variability could be used to evaluate soil moisture predictions from more 
sophisticated, complex soil moisture models. As expected, evergreen systems showed reduced 
seasonal and interannual variations. This behavior was true for precipitation, but our results showed 
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that 50% or more of the NDVI variance in May can be attributed to temperature variability (Figure 9c) 
for some regions in the boreal forests. Thirty to sixty percent (30-60%) of the NDVI variance for the 
Sahel-Sahara border of northern Africa and in Southern Africa can be explained by precipitation 
variability during their respective rainy seasons. During their respective diy seasons NDVI variance 
was low in general. Temperate grasslands showed strong climate (precipitation) controls during the 
growing season. Here again NDVI variance tended to be low outside the growing season so low 
explained variance by climate was expected for these periods. The climate driven variability in NDVI 
shown in this study demonstrates the ability of the GIMMS NDVI3g data record to provide realistic 
representation of interannual variability in solar radiation absorbed by vegetation. 
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